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Before the availability of digital photography resulting from the invention of charged couple devices in 1969, the measurement 
of plant architecture was a manual process either on the plant itself or on traditional photographs. The introduction of cheap digital 
imaging devices for the consumer market enabled the wide use of digital images to capture the shape of plant networks such as 
roots, tree crowns, or leaf venation. Plant networks contain geometric traits that can establish links to genetic or physiological 
characteristics, support plant breeding efforts, drive evolutionary studies, or serve as input to plant growth simulations. Typically, 
traits are encoded in shape descriptors that are computed from imaging data. Skeletons are one class of shape descriptors that are 
used to describe the hierarchies and extent of branching and looping plant networks. While the mathematical understanding of 
skeletons is well developed, their application within the plant sciences remains challenging because the quality of the measure- 
ment depends partly on the interpretation of the skeleton. This article is meant to bridge the skeletonization literature in the plant 
sciences and related technical fields by discussing best practices for deriving diameters and approximating branching hierarchies 
in a plant network. 
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The advent of digital imaging accelerated the development 
of computational methods to digitally reconstruct the shape of 
plant networks from imaging data. From reconstructions, de- 
tailed measurements such as network width or diameters of 
network links can be derived for a plant network. The impor- 
tance of reconstructions and measurements (Godin, 2000; 
Spalding and Miller, 2013) is founded in the many known or 
hypothesized links between plant network properties and ac- 
tual functions of plants (Brodribb and Feild, 2000; Enquist, 
2002; Lynch, 2011; Lynch and Brown, 2012). In particular, 
roots system architecture (Lynch, 1995; de Dorlodot et al., 
2007), tree crown architecture (Delagrange et al., 2004), or 
leaf vein structures (Roth-Nebelsick et al., 2001; Mason et al., 
2013; Price et al., 2013) are examples of plant networks. Mea- 
surements of such networks are analyzed in the context of 
genotype to phenotype mapping (El-Soda et al., 2014) as well 
as breeding efforts (Holland et al., 2003; Finkel, 2009), evolu- 
tionary studies of plant species (Green and Hickey, 2005; 
Fernandez et al., 2010; Feild et al., 2011), or plant growth 
simulation (Vansteenkiste et al., 2014). 

Images encode the shape of plant networks such as roots, 
tree crowns, or leaf venation. Most commonly, images are 
digitally represented in a two- or three-dimensional raster of 
pixels or of voxels, respectively. Some imaging techniques 
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produce point clouds of 3D coordinates that often employ no 
neighborhood information between the points. Recently, the 
combination of imaging data with skeletonization algorithms 
became popular as a tool to extract the hierarchies and geo- 
metric measurements of plant networks. Skeletons consist of 
a collection of (approximated) curves being mostly centered 
within the plant network of interest (Fig. 1). In this article, I 
discuss the fundamental characteristics of imaging data in the 
context of studying plant networks and give an overview of 
existing skeletonization concepts. Within the skeletonization 
overview, I highlight concepts that were applied to plant net- 
works. In the remainder of the article, I introduce two basic 
skeleton concepts in detail and demonstrate their applicability 
to plant networks on selected examples. 

THE BOUNDARY CONCEPT AND OCCLUSION IN 
IMAGING DATA 

Before applying skeletonization algorithms, it is important 
to understand the intrinsic constraints of imaging data. Imag- 
ing technologies capture data about an object surface and/or 
interior within the specification of an imaging sensor. For ex- 
ample, a digital photo camera collects the light reflected from 
surfaces that face the photo sensor. At each location of the 
photo sensor, a value for the received wavelength is stored as 
an electrical potential or as the number of photons counted 
within a certain time frame. The stored wavelength corre- 
sponds to colors that we recognize in the pixels of the image. 
In the case of a plant network, this physical process of captur- 
ing light is a projection of the plant network onto a plane. The 
result is a 2D image that encodes the form of a plant network. 
For leaf venation networks, this projection does not cause any 
self-occlusion (Fig. 2). However, the projection onto the 2D 
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Fig. 1. 3D point cloud data obtained with a terrestrial laser scanner and the corresponding curve approximation skeleton of an orchard tree. The skeleton 
was computed with SkelTre (Bucksch et al., 2010). 



plane comes at a cost for other plant networks such as heavily 
branched tree crowns or dense roots. In both cases, branches 
will occlude parts of the tree crown (Fig. 2) or rootstock. As a 
result, the digital photograph misses network information at 
all self-occluded locations. Similar to recording reflections of 
light, it is also possible to record absorption rates of wave- 
lengths that penetrate the plant network. 

3D imaging technologies improve upon the occlusion prob- 
lem. In general, 3D data are obtained from the reflection of 
waves from the plant surface or from waves penetrating the 
plant network such that the absorption rate varies through dif- 
ferent tissues. Technologies relying on the reflection of light 




leaf venation tree crown 

Fig. 2. Plant networks with and without occlusion. (A) Image of a 
cleared leaf with its venation network visible. The venation network is not 
occluding itself through the projection onto the image plane. Image taken 
from http://clearedleavesdb.org (Das et al., 2014). (B) Image of a tree 
crown. The branches occlude each other, such that the image does not con- 
tain the information of the complete branches. 



from a surface include structured light methods that project a 
known pattern onto the plant network and calculate 3D infor- 
mation from the distortion of the projected pattern (Salvi 
et al., 2004). Another technology is terrestrial laser scanning 
(Vosselman and Maas, 2010) that relies on the reflection of 
laser beams from the surfaces. The round- trip time of the laser 
from the sensor to surface determines the distance between 
sensor and surface because of the constant speed of light. In 
practice, occlusion effects still occur as missing data, but can 
be minimized by using registration techniques to combine dif- 
ferent views of the plant network (Bucksch and Koshelham, 
2013). The benefit of the obtained 3D information is that dis- 
tribution of the network in 3D space allows the use of prior 
knowledge, such as assuming allometric relations in a branch- 
ing structure (Aiteanu and Klein, 2014), to recover an approx- 
imation of the network. Ideally, the occlusion problem is 
resolvable with imaging technologies that penetrate the plant 
network, such as X-ray (Mairhofer et al., 2013) or magnetic 
resonance imaging (Schulz et al., 2013), that record changing 
absorption rates in different plant tissues. As a result, not only 
the outer surface, but also the internal architecture of the plant 
is resolvable. The important commonality between all these 
imaging technologies is that the plant network is represented 
as a boundary between itself and the free space surrounding 
the object in 2D and 3D. Or, more intuitively spoken, the 
boundary of a 3D plant network separates its enclosed volume 
from its "outside" environment. Data representations that just 
encode the plant network and the "outside" are called binary, 
because they assign a 1 to every plant network location and a 
0 to every background location. In most application scenarios, 
we can safely assume the boundary of a plant network to be 
smooth, i.e., the boundary has no sharp corners or breaklines. 

Once we understand the concept of a boundary, we can ob- 
serve that skeletons are derived from this interface between the 
plant network and environment. Thus, a skeleton representation 
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is computed on the basis of the boundary and, therefore, a 
skeleton represents boundary properties rather than the hierar- 
chies and extensions of biological interest. 



SKELETON CONCEPTS AND PLANT APPLICATIONS 

Current literature reviews a plethora of skeleton concepts to 
be computed from 2D binary images (Lam et al., 1992), 3D im- 
ages of voxels (Sobiecki et al., 2014), and point clouds (Cornea 
and Min, 2007; Bucksch, 2011). Two basic skeleton concepts 
used today are the medial axis (Lee, 1982) and the Reeb graph 
(Biasotti et al., 2008) shown in Fig. 3. The computation of the 
medial axis can be performed on 2D binary images (Serra, 
1986), 3D voxels (Palagyi, 2002), and point clouds (Amenta 
et al., 2001). A noteworthy extension is the Hamilton- Jacobi 
skeleton (Siddiqi et al., 2002), which enhances the medial axis 
with a flux equation to achieve more robustness to noise. In re- 
cent years, the field of computer graphics developed several 
new skeleton concepts that make use of prior knowledge of the 
plant network. Rotational symmetry axis (ROSA) assumes ro- 
tational symmetry around the skeletal curves to enhance robust- 
ness against missing data from point clouds (Tagliasacchi et al., 
2009). The assumption of a tree structure to extract a skeleton 
from point cloud data was introduced in Livny et al. (2010) for 
trees and has been applied to plant networks of grape clusters 



(Scholer and Steinhage, 2012; Steinhage et al., 2012). A very 
recent concept is the L r medial skeleton (Huang et al., 2013), 
which shows promising results on point clouds of plants. 

The use of skeletonization algorithms to obtain quantita- 
tive measurements of plant networks was first applied on tree 
crowns to estimate the radiative properties by simulating the 
foliage and fine branches from laser data (Cote et al., 2009) 
and to derive branch hierarchy information along with length and 
diameter distributions (Bucksch and Fleck, 2011; Delagrange 
et al., 2014). Leaf veins were quantified with skeletons to 
obtain estimates of areole measures and hierarchies (Price 
et al., 2011, 2013, 2014; Price, 2012). A recent article intro- 
duced a framework to identify legumes on the basis of leaf 
vein features (Larese et al., 2014). The branching structure of 
roots is another plant network that is often assessed with 
skeletons to obtain measurements of phenotypic traits (Iyer- 
Pascuzzi et al., 2010; Clark et al., 201 1 ; Galkovskyi et al., 2012; 
Pound et al., 2013) to be related to regions of the genome 
(Topp et al., 2013). Pound et al. (2013) introduced a skeletoni- 
zation concept for plant phenotyping that uses a probability 
field computed from the gray values of a 2D image to grow a 
connected skeleton of young and sparse roots. The growth is 
initiated at a user-defined location to the root tips. Free soft- 
ware exists that allows manual enforcement of skeleton cen- 
teredness and correction of miscomputed branching hierarchy. 
The references for tree crowns are examples of successful 




Fig. 3. Comparison of the medial axis and the Reeb graph for a simple rectangle. The medial axis (MA) was extracted using a morphological thinning 
algorithm. Examples are shown of corresponding circles whose center (c) build the medial axis and touch the boundary of the plant network M in two 
points. The Reeb graph (RG) was computed using the height function F to extract the contours (green line) at height levels i x and i 2 from the rectangular 
plant network M at intervals of s. The center of the contracted contour is marked with a cross. The source code for all examples is available on 
https://github.com/abucksch. 
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undistorted boundary 



boundary distorted by 2 pixels 



Fig. 4. Sensitivity of the medial axis to boundary distortions. Two pixels were added to the boundary, which results in two extra branches that do not 
belong to the plant network hierarchy. 



applications developed on the basis of a Reeb graph. Except 
for Pound et al. (2013), the chosen examples for leaves and 
roots apply to the medial axis. 

Medial axis — The medial axis (Fig. 3) was introduced by 
Harry Blum (Blum, 1967) and was inspired by line-like 
structures created by the process of extinguishing fires. Simply 



imagine a patch of dry grassland, the whole border of which 
is ignited at one moment in time. If the fire propagates in- 
ward at a constant speed, then the fire fronts meet at the loca- 
tion of the medial axis. Although there is no generally 
agreed-upon definition of the medial axis (Attali et al., 
2009), we can use here a simplified definition that can be 
visually followed in Fig. 3 using the same variables. The medial 



undistorted boundary 



boundary distorted by 2 pixels 



Fig. 5. Reeb graph of a simple branching structure extracted with the height function. Vertices are marked in red. Edges are marked in green if the 
distance between two vertices is larger than two pixels. At these center points, diameter measurements can be derived without interpolation of a curve be- 
tween two vertices. The level sets were extracted at intervals of one pixel. The robustness of the Reeb graph to small disturbances of the boundary is dem- 
onstrated in the example in Fig. 3. 
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axis MA of a plant network M with a smooth boundary is the 
set of points inside the root network that consists of centers 
(c;) of circles in two dimensions or spheres in three dimen- 
sions, whose boundary touches in at least two points (t iyj and 
t if j +1 ) the boundary of the plant network. Here, touching is 
used as a synonym for the circle/sphere being tangent to the 
boundary of M. 

The definition of the medial axis implies that each local 
change of the boundary changes the medial axis. Such changes 
correspond to locally occurring/vanishing maxima, minima, or 
saddle points on the boundary. Consequently, additional branches 
occur in the medial axis that may not represent the actual net- 
work. Reasons for additional branches include noise in the im- 
aging data and tissue growth that results in a locally bulging 
boundary (Fig. 4). As a result, direct biological interpretation 
from the medial axis skeleton is problematic. A simple but ef- 
ficient correction that makes use of a property of the medial axis 
is to remove branches that are shorter than the known radius of 
a circle/sphere at a branching point plus a user-defined thresh- 
old. In some applications, a threshold might be hard to define, 
such that methods to prune the medial axis have to be applied 
(e.g., Bai et al., 2007). 

One problem arising in 3D applications of the medial axis is 
the occurrence of planes. For many plant applications, this 
problem is irrelevant, because plant organs are often cylindrical 
or planar and therefore result in curves. However, there are so- 
lutions to reduce these planes to curves (Palagyi and Kuba, 
1998; Gorte and Pfeifer, 2004). 



Reeb graph — The Reeb graph (Fig. 3) is a skeletonization 
concept introduced in 1946 (Reeb, 1946) and used to describe 
the change of an object shape in relation to a function ex- 
panded over the object. A simple way to imagine a Reeb graph 
is to think of a banana sliced — lengthwise — by a knife. The path 
of the knife along the banana can be thought of as the function 
that is expanded over the banana. The centers of the cutting 
area through the banana are the points on the Reeb graph that 
get connected to the center of the next cutting area in direction 
of the path. The following formal definition can be graphically 
followed in Fig. 3. Given a plant network M with smooth 
boundary and a piecewise linear function F, the level set at a 
value i is defined as the set of points in M with function value 
equal to i. As a practical example, if F is the height function 
whose function values are the values along the y-axis of a Car- 
tesian coordinate system, then the level sets are simply the 
intersections when slicing through the object at intervals of 
size s orthogonal to the y-axis. To obtain the Reeb graph, we 
analyze the change of the connected components of the level 
sets of F, called contours. The Reeb graph RG of F is obtained 
by contracting the contours of F to points. For plants, this 
point can be chosen as the center of the contour, because most 
contours are convex. The center points of the contours assem- 
ble the vertices of the Reeb graph, which are connected via 
edges to the "closest" vertices obtained at the next function 
value i+s. For easier understanding, the given definition in- 
cludes the embedding of the graph into the plant network. The 
original literature on Reeb graphs uses an abstract definition, 
which only defines the merge and split of contours free of co- 
ordinates. An example of the introduced construction is shown 
in Fig. 5 with the same binary images used in Fig. 4. Addition- 
ally, Fig. 5 demonstrates the robustness of the Reeb graph to 
small changes on the boundary. 



The contours obtained by using height function are horizon- 
tal slices through the plant network. Problems can arise because 
the height function is not rotation invariant. To understand the 
potential problems, consider the Reeb graph of a vertical trunk 
of a plant network with a child branch emerging at a 90° angle 
from a parent branch (Fig. 6). At the level of the child branch, 
the contours include both the child and parent branch. Hence, 
only one contour for parent and child branch is extracted. As 
a result, consecutive center points appear to be "far from each 
other" or "shifted." 

A different choice of the function F can result in different 
Reeb graphs that can overcome the problems of the height func- 
tion. Several solutions have been suggested to choose a suitable 
function F. Known solutions include the choice of different 
global functions (Biasotti et al., 2000) or determining a suited 
function locally (Biasotti et al., 2008; Bucksch et al., 2010), as 
shown in Fig. 7. 

APPROXIMATING THE BRANCHING HIERARCHY OF A 
PLANT FROM REAL DATA 

In Fig. 8, the medial axis and its segmentation of a cowpea 
root system was computed from a 2D image. The segmenta- 
tion allows the detection of branches that do not represent the 
actual plant architecture. The segmentation was achieved by 




Fig. 6. Problems arising from the choice of the contour-extracting 
function F. The height function is not rotation invariant, which results in 
a local shift of the Reeb graph at the location where a branch is emerging 
at a 90° angle. 



http://www.bioone.org/loi/apps 



5 of 10 



Applications in Plant Sciences 2014 2(8): 1400005 
doi:10.3732/apps. 1400005 



Bucksch — A practical introduction to skeletons for the plant sciences 




J 



Fig. 7. Solution to rotation invariance. (A) A point cloud of a tree obtained with a terrestrial laser scanner. (B) A computed SkelTre skeleton that con- 
tains the Reeb graph. The computation changes locally the function to extract the contours. (C) The extracted contours of the tree on the left colored ran- 
domly. The highlighted red boxes show a noncentered skeleton part and the extracted contours. Image taken from Bucksch (201 1). 



tracing the shortest path from the tip of the medial axis closest 
to the top row of the image to every root tip that is represented 
in the medial axis. The tracing resolves approximately for the 
loops resulting from the projection of the three-dimensional 
arrangement of roots to the image plane. Note that this strat- 
egy only applies to sparse networks, where branches do not 
occlude each other. 



Figure 9 shows a 3D skeleton of a tree crown. Here, the 
skeleton is represented as a graph embedded into the original 
point cloud of a tree. All vertices with one incident edge that 
connect to vertices with three or more incident edges are re- 
moved from the originally extracted graph. The resulting skel- 
eton approximates the actual branching structure of the tree 
crown (Bucksch, 2011). 




Fig. 8. Detecting medial axis parts for removal. (A) A digital image of a cowpea root. (B) The medial axis derived from the image. (C) The segmented 
skeleton with random colors for each detected branch. Smaller branches are now removable by using a threshold and the known radius of the medial circle 
at the tip of each branch. The green box indicates an example area where additional branches can be removed. 
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SKELETON CENTEREDNESS USED TO DERIVE 
DIAMETERS IN PLANT APPLICATIONS 

For plant networks, noncentered skeleton parts occur at loca- 
tions where a child branch emerges from a parent branch (Fig. 7 
and Fig. 10). Noncenteredness is implicit because a skeleton is 
a one-dimensional representation of a plant network that ex- 
pands in two or three dimensions. This expansion is what we 
visually recognize as thickness and can be seen as the result of 
the expansion process of a growing one-dimensional curve 
(Bucksch et al., 2014). Therefore, every connected skeleton of 
a branching plant network travels from a branching point 
through a region of "thickness" before it reaches the location 
where the branch emerges from the boundary (Fig. 10). As a 
result, the skeleton part between the branching point and the 
point of branch emergence on the boundary is not centered. In 
the case of the medial axis, we can simply subtract the radius of 
the medial circle at the branching point from the emerging skel- 
etal branch. In the case of the Reeb graph, we may identify jumps 
in the size of the contour on suited data. As a result, we obtain an 
approximation of the point of branch emergence on the bound- 
ary. Computing the point of branch emergence corresponds 
closely to measurements that can be obtained manually on the 
real plant, e.g., with a measuring tape, when measuring along the 
surface of the branch from the start of a branch on the trunk sur- 
face to the tip. For the Reeb graph, we can detect such regions as 
"long edges." Note that vertices are favorable locations to derive 




(B) Vertices with only one incident edge and an adjacent vertex with three 



measurements from the plant network because the graph can 
maintain centeredness for its vertices while the edges are not 
necessarily centered. 



A QUANTITATIVE EXAMPLE OF SKELETON-DRIVEN 
TRAIT COMPUTATION 

In the following, I apply the introduced skeleton concepts to 
images from a freely available rice root data set used in Iyer- 
Pascuzzi et al. (2010). The data set consists of images that show 
the root from several side views. I selected one image per indi- 
vidual root of the Caipo and Bala genotypes. Each image was 
selected such that occlusion of roots is visually minimized, 
which corresponds to the view of maximal root width. Overall, 
I used nine images of the Caipo data set and eight images of the 
Bala data, shown in Fig. 11. Both data sets were imaged with 
the same imaging setup such that results measured in pixels are 
comparable. The origin of the image is set to be the upper left 
corner. The absolute value of the negative y-direction corre- 
sponds to the row numbers. 

The three example traits — maximum tip diameter, root 
width, and rooting depth — were computed for the selected Bala 
and Caipo images. Caipo and Bala were chosen because they 
visually contrast with each other in these traits (Fig. 11). A tip on 
the medial axis is defined as a pixel with only one neighbor and 
the diameter is known from the corresponding medial circle. To 
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Fig. 10. Noncentered locations of the medial axis within the green 
medial circle. To be a connected skeleton, the medial axis passes through 
an area of the object interior before reaching the point of branch emergence 
marked with a cross. 



estimate root width, I used the difference between maximal and 
minimal column values occupied by the skeleton. Rooting 
depth corresponds to the maximal row number reached by one 
of the tips. I also applied noise filtering with a threshold of two 
pixels to eliminate branches resulting from noise as described 
above. 

Caipo shows a maximum tip diameter of 3.9 pixels with a 
standard deviation (SD) of 0.57. Bala resulted in a maximum 
tip diameter of 4.4 pixels (SD = 0.50). A Mest performed on the 
two measurement sets resulted in a P = 0.13 and suggests that 
the maximum tip diameter differentiates both genotypes. Fur- 
thermore, I extracted estimates of the two known traits — root 
width and rooting depth. Root width resulted in an average of 
607.1 1 pixels (SD = 224.90) for Caipo and an average of 504.38 
pixels (SD = 66.14) for Bala. Although the r-test resulted only 
in P = 0.23 for the small set of test images, the genotypes dif- 
fered strongly in their standard deviation. The strongest differ- 
entiation (P < 0.0001) between the two genotypes was achieved 
with the rooting depth. Rooting depth of Caipo was 1317.80 
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pixels (SD = 151.35) on average. In comparison, the rooting 
depth of Bala was about half as deep as the rooting depth of 
Caipo, resulting in an average of 625.50 pixels (SD = 181.62). 
I used an early version of the Digital Imaging of Root Traits 
(DIRT) platform (www.dirt.biology.gatech.edu) to compute 
the traits from the medial axis. The computation results are in- 
cluded in Appendix S 1 . 



DISCUSSION 

The advancement of image processing algorithms has re- 
sulted in a wide variety of software tools throughout biology 
(Eliceiri et al., 2012) and the plant sciences in particular (Lobet 
et al., 2013). In addition to end-user software, open-source li- 
braries for scripting languages such as MATLAB or Python 
(Coelho, 2013; van der Walt et al., 2014) offer access to skele- 
tonization methods. In particular, Python and its scientific li- 
braries have become widely used in interdisciplinary research 
and education (Deutsch, 2014). Software libraries enable re- 
searchers to adapt skeletonization methods to their specific ap- 
plications. Specific application needs can often be solved via 
open community platforms such as MATLAB Central, stack 
overflow, or stack exchange that offer access to source code for 
easy tasks, such as implementations of 3D skeletonizations in 
MATLAB (Schena, 2004) or finding branches in a skeleton 
(stackoverflow.com, 2013). From a user perspective, I intro- 
duced two skeletonization basic concepts and provided insight 
into mechanisms that potentially influence the statistical analy- 
sis of data computed with available software packages or code 
snippets collected on the Internet. However, it remains the re- 
sponsibility of the researcher to judge and document if mea- 
surement errors induced by imperfect data or the skeletonization 
concept have significant influence on the achieved results. 

CONCLUSIONS 

Skeletons are powerful descriptors for analyzing plant net- 
works, if applied with care. In this article I pointed out two 
difficulties: 

1 . Spurious branches have to be removed from an extracted 
skeletal structure to obtain an improved approximation of the 
actual plant network. 

2. Diameter measurements of the plant network of interest 
cannot be derived at all locations of the skeleton. 

The choice for a particular skeletonization concept is depen- 
dent on the application, but knowledge of the basic underlying 
concepts will make it easier to adapt source code to plant ap- 
plications. These adaptions result in improved data quality 
computed from skeletons. As a consequence of improved data 
quality, the statistical analysis of the naturally large geometric 
variations within plant networks is facilitated. 

Achieving a one-to-one correspondence between the actual 
plant network and a skeleton extracted from imaging data are 
only possible on sparse plant networks. The limiting factor is not 
the skeleton concept, but the imaging technology and its ability 
to resolve the complete plant network. Therefore, preparing near- 
ideal input data will involve manual work to correct for the phys- 
ical limits imposed by certain imaging technologies. The simple 
examples shown here should help end-users to judge the results 
computed with the chosen software or encourage users to contact 
the author of the software to interpret the results correctly. 
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Bala 



Caipo 




Fig. 11. The Bala and Caipo images. Maximum tip diameter, root width, and rooting depth were extracted to differentiate both genotypes from 
each other. 



Although I used several thresholds in this article, I do not 
propose them as an optimal algorithmic solution, but as a sim- 
ple tool to start experiencing and understanding skeleton con- 
cepts. Appendix S2 contains annotated references to freely 
available data, source/pseudo codes, and executable software 
referenced in this article. 
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